

package MGD;

use strict;
use base qw(Exporter);

use List::Util qw(reduce);
use Math::Matrix;

our @EXPORT_OK = qw(
                    ComputeMGD
                );

our %EXPORT_TAGS = (
    all => [ @EXPORT_OK ],
);

use constant PI => 4 * atan2(1,1);

#
#Compute the mean value for a vector
#Parameters:
#    vector: reference to an input vector (list)
#Output:
#    the mean value for that vector
sub mean_vector {
    my $vector = shift;

    return 0 if (@$vector == 0);

    return (reduce {$a + $b} @$vector) / scalar @$vector;
}


#Compute the mean value for a matrix
#Parameters:
#    matrix: an m * n Math::Matrix,
#            m is the number of samples, n is the number of features
#Output:
#    a 1 * n one diminsion Math::Matrix
#        each value is the the mean for each feature in the matrix
sub mean_matrix {
    my $matrix = shift;

    my ($m, $n) = $matrix->size();

    return () if ($m == 0 || $n == 0);

    my $t_matrix = ~$matrix;

    return Math::Matrix->new([map { mean_vector($t_matrix->[$_])} 0 .. $n-1]);
}


#Compute the probility
#Parameters:
#    matrix: an m * n array, the parameter must be a reference
#            m is the number of samples, n is the number of features
#Output:
#    reference of an m * (n+2) array,
#        the first column is the probility,
#        the second column is the offset divided by the standard deviation,
#               if it is larger than 4, likely to be anomaly ( < 0.1% )
#        others are the original array
#
sub ComputeMGD($) {
    my $array_ref = shift;
    my $x = Math::Matrix->new(@$array_ref);
#print print "Matrix:\n", $x;

    my $mu = mean_matrix($x);
#print print "Mu:\n", $mu;

    my ($m, $n) = $x->size;

    my $sigma = Math::Matrix->new(map { [map {0} 1 .. $n] } 1 .. $n ); # Zero(n * n)

    for my $i (0 .. $m-1) {
        my $xi = Math::Matrix->new($x->[$i]);
        $sigma += ~($xi - $mu) * ($xi - $mu);
    }
    $sigma = $sigma->multiply_scalar(1/$m);
#print "Sigma:\n", $sigma;

    my $sigma_inv = $sigma->invert() || $sigma->pinvert();
#$sigma_inv->print("SigmaInv:\n");

    #the standard deviation
    my $sd = sqrt($sigma->determinant());
#print "sd: $sd\n";
    my $factor = $sd == 0 ? 1 : 1 / ( ((2 * PI) ** ($n / 2)) * $sd );

    #compute p(x) for each x
    my @probs = ();
    for my $i (0 .. $m-1) {
        my $xi = Math::Matrix->new($x->[$i]);
        my $xi_mu = $xi - $mu;

        my $dist = $xi_mu->absolute(); #offset from mu

        my $product = ($xi_mu * $sigma_inv)->dot_product($xi_mu);


        push @probs, [ $factor * exp(-0.5 * $product), 
                       $sd == 0 ? 0 : $dist/$sd ];
    }

    return \@probs;
}


1;
